$$ A \left({{a}^\dagger} + {a}\right) + B {\sigma_x} + C \left({{b}^\dagger} + {b}\right) - \Delta_{d} {{a}^\dagger} {a} - \Delta_{p} {{b}^\dagger} {b} + \Delta_{s} {\sigma_z} + g \left({\sigma_-} {{a}^\dagger} + {\sigma_+} {a}\right) + \lambda \left({{b}^\dagger} + {b}\right) {\sigma_x} $$

In [1]:
#import functions
%pylab inline

# from MyUnits import *
from MyFunctions import *
from qutip import *

# from MyQubit import *
# import mpld3
import multiprocessing as mp
import itertools


Populating the interactive namespace from numpy and matplotlib

In [2]:
import scipy.constants as sc

In [3]:
import time

In [4]:
def calc_spectrum_3(N, M, Delta_C,Delta_NR,
                    Delta_Q, g, lamb, 
                    A_C, A_Q, A_NR,
                    kappa_n, gamma_rel,
                    gamma_dep,n_th_a):
    """
    Calculate Qubit density Matrix rho and trace of: rho_sz, rho_n, rho_x, rho_a : Based on work of : Suri, B. et al. Observation of Autler–Townes effect in a dispersively dressed Jaynes–Cummings system. New J. Phys. 15, 125007 (2013).
    look also calc_spectrum function
    N: number states cavity
    t_list: time - not used
    Delta_d_tilde: cavity tone frequency difference
    Delta_s_tilde: qubit tone frequency difference
    chi: dispersive coupling
    Gamma_d: cavity tone power (GHz)
    Gamma_s: qubit tone power (GHz)
    """


    # cavity operators
    a = tensor(destroy(N), qeye(M), qeye(2))
    n_C = a.dag() * a
    x_C = a + a.dag()

    # NR operators
    b = tensor(qeye(N), destroy(M), qeye(2))
    n_NR = b.dag() * b
    x_NR = b + b.dag()

    # qubit operators
    sm = tensor(qeye(N), qeye(M), destroy(2))
    sz = tensor(qeye(N), qeye(M), sigmaz())
    sx = tensor(qeye(N), qeye(M), sigmax())
    nq = sm.dag() * sm
    xq = sm + sm.dag()

    # Unitary
    I = tensor(qeye(N), qeye(M), qeye(2))


    # Hamiltonian

    H_C = Delta_C * (a.dag() * a + I / 2.0)  # Cavity

    H_NR = Delta_NR * (b.dag() *b + I / 2.0 ) # Nanoresonator

    H_Q = (Delta_Q / 2.0) * sz  # qubit

    H_QtoC = g * (sm * a.dag() +  sm.dag()*a )  # interaction Cavity Qubit

    H_QtoNR = lamb * (b.dag() + b)*sx # interaction NR Qubit

    H_A = A_C * (a.dag() + a) # tone cavity

    H_B = A_Q * (sm.dag() + sm) # tone qubit

    H_C = A_NR * (b.dag() + b) # tone NR



    # # extra term Unitary tranformation and RWA
    # He = nc * chi ** 2 / (4 * Delta) * (w_q / Delta - I)
    # Hd = (Gamma_d / 2.0) * (a + a.dag())  # Drive cavity
    # Hs = (Gamma_s / 2.0) * (sm + sm.dag())  # Drive qubit
    # HKerr = zeta_1 * ((a.dag() * a) ** 2) * sz + zeta_2 * (a.dag() * a) ** 2  # non linear Kerr terms

    H = H_C + H_NR + H_Q + H_QtoC + H_QtoNR + H_A + H_B + H_C  # total Hamiltonian

    # collapse operators

    c_op_list = []

    rate = kappa_n * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * a)

    rate = kappa_n * n_th_a
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * a.dag())

    rate = gamma_rel * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * sm)

    rate = gamma_rel * (n_th_a)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * sm.dag())

    rate = (gamma_dep / 2) * (1 + n_th_a)
    if rate > 0.0:
        c_op_list.append(sqrt(rate) * sz)

    # Calculation stedy state Master equation
    rho = steadystate(H, c_op_list)

    # variables to be returned

    rho_sz = rho * sz
#     rho_n = rho * (a.dag() * a + I / 2)
#     rho_x = rho * (a.dag() + a)
    rho_a = rho * a
    rho_b = rho * b

    return rho_sz.tr(), rho_a.tr(), rho_b.tr()

In [5]:
N = 5
M = 5

Delta_C = 1
Delta_NR = 1
Delta_Q =  1

g = 0.119
lamb = 0.1

A_C = 0.000308
A_Q = 0.00026
A_NR = 0.000001

kappa_n = 0.00025
gamma_rel = 0.00023
gamma_dep = 0.00114
n_th_a = 0.3

In [6]:
def Delta(f_r,f_d):
    """
    f_r resonance freq
    f_d drive freq
    """
    return f_r-f_d

def Gamma(P,Q_c,Q_l,kappa,w):
    return sqrt(Q_c*kappa*1e9*dBmtoW(P)/(2*Q_l*sc.hbar*w*1e9))/1e9

In [ ]:


In [7]:
# a1, b1 = zip(*itertools.product(Delta_C, Delta_Q))
# pool = mp.Pool(processes=12)
# results2 = [pool.apply_async(calc_spectrum_3,
#                             (N,
#                              M,
#                              C,
#                              Delta_NR,
#                              Q,
#                              g,
#                              lamb,
#                              A_C,
#                              A_Q,
#                              A_NR,
#                              kappa_n,
#                              gamma_rel,
#                              gamma_dep,
#                              n_th_a)) for C,Q in zip(a1,b1)]

# results = [p.get() for p in results2]
# results.sort()

Device Data


In [8]:
Ej = 11.55 # Maximum Josephson Energy

Ec = 0.22 # Capacitive Energy

w_ge_max = sqrt(8*Ec*Ej)-Ec
print('Qubit max frequency: ',w_ge_max, 'GHz')


f = w_ge_max *1e9 # GHz
T = 50e-3 # K

n_th_a = 1/(exp(sc.h*f/(sc.k*T)-1))

print('<n> thermal: ',n_th_a)


Qubit max frequency:  4.28865833702 GHz
<n> thermal:  0.0443136288802

In [10]:
Q_c = 20000
Q_l = Q_c


N = 5
M = 5


g = 0.1
lamb = 0.00011



kappa_n = 0.0002
gamma_rel = 0.0002
gamma_dep = 0.001


# Cavity
w_c = 5 * (2 * pi)
w_d = w_c + 0.006283185307175643 #(2 * pi) * 4.995
P_C = -120# linspace(-60,-150,20)

# Nanoresonator
w_NR = 3 * (2 * pi)
w_p = (2 * pi) * 2.9996
P_NR = -1000


# Qubit
w_q =  4.5 * (2 * pi)
w_s = (2 * pi) * linspace(4.2,4.8,300)#4.501
P_Q = linspace(-50,-150,50)

In [11]:
Delta_C = Delta(w_c, w_d) 

Delta_NR = Delta(w_NR, w_p)

Delta_Q = Delta(w_q, w_s)
# print(Delta_C,Delta_NR,Delta_Q)

In [12]:
A_C = Gamma(P_C,Q_c,Q_l,kappa_n,w_c)

A_Q = Gamma(P_Q,Q_c,Q_l,kappa_n,w_c)

A_NR = Gamma(P_NR,Q_c,Q_l,kappa_n,w_c)
# print(r'A_C = ', A_C, '\n\nA_Q = ', A_Q, '\n\nA_NR = ' , A_NR)

In [13]:
a1, b1 = zip(*itertools.product(A_Q, Delta_Q))

Simulation


In [ ]:
pool = mp.Pool(processes=12)
t0 = time.time()
results1 = [pool.apply_async(calc_spectrum_3,
                            (N,
                             M,
                             Delta_C,
                             Delta_NR,
                             y,
                             g,
                             lamb,
                             A_C,
                             x,   
                             A_NR,
                             kappa_n,
                             gamma_rel,
                             gamma_dep,
                             n_th_a)) for x , y in zip(a1,b1)]

t1 = time.time()
print("elapsed =", (t1-t0))
# results2 = [pool.apply_async(calc_spectrum_3,
#                             (N,
#                              M,
#                              Delta_C,
#                              Delta_NR,
#                              Delta_Q,
#                              g,
#                              lamb,
#                              A_C,
#                              A_Q   
#                              A_NR,
#                              kappa_n,
#                              gamma_rel,
#                              gamma_dep,
#                              n_th_a)) for x , y in zip(a1,b1)]

results = [p.get() for p in results1]

In [ ]:
results = [p.get() for p in results1]

In [ ]:
shape(results)

In [ ]:
Tr_sz = reshape(array(results)[:,0],(-1,len(Delta_Q+1)))
Tr_a = reshape(array(results)[:,1],(-1,len(Delta_Q+1)))
Tr_b = reshape(array(results)[:,2],(-1,len(Delta_Q+1)))

In [ ]:
fig, ax = subplots(2,2, figsize=(16,10))

im = ax[0,0].pcolor(  w_s/2/pi,P_Q,abs(Tr_sz))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0,0])
# ax[0,0].set_xlim(4.27,4.39)
# ax[0,0].set_ylim(P_i,P_f)
ax[0,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0,0].set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)


im2 = ax[0,1].pcolor(w_s/2/pi,P_Q,abs(Tr_a))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[0,1])
# ax[0,1].set_ylim(P_i,P_f)
# ax[0,1].set_xlim(4.27,4.39)
ax[0,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[0,1].set_title(r'$Tr[\rho {a}]$',fontsize=20)

im3 = ax[1,0].pcolor(w_s/2/pi,P_Q,abs(Tr_b))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,0])
# ax[1,0].set_ylim(P_i,P_f)
# ax[1,0].set_xlim(4.27,4.39)
ax[1,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,0].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,0].set_title(r'$Tr[\rho (b)]$',fontsize=20)

In [ ]: